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Abstract 

■ Dissipative Particle Dynamics (DPD) is becoming a popular particle based method to study 

flow through microchannels due to the ease with which the presence of biological cells or DNA 
chains can be modeled. Many Lab-On-Chip (LOG) devices require the ability to manipulate the 
^ ■ transport of cells or DNA chains in the fluid flow. MicroChannel surfaces coated with combinations 

. of hydrophilic and hydrophobic materials have been found useful for this purpose. In this work, 

we have numerically studied the hydrodynamics of a steady nonuniform developing flow between 
two inflnite parallel plates with hydrophilic and hydrophobic surfaces using DPD for the flrst time. 
The hydrophobic and hydrophilic surfaces were modeled using partial-slip and no-slip boundary 
conditions respectively in the simulations. We also propose a new method to model the inflow and 
outflow boundaries for the DPD simulations. The simulation results of the developing flow match 
analytical solutions from continuum theory for no-slip and partial-slip surfaces to good accord. 

The entrance region constitutes a considerable fraction of the channel length in miniaturized 
devices. Thus it is desirable for the length of the developing region to be short as most microfluidic 
devices such as cell or DNA separators and mixers are designed for the developed flow field. We 
studied the effect of a hydrophilic strip near the inlet of a microchannel on the effective developing 
length. We find that the presence of the hydrophobic strip significantly reduces the developing 
length. 
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I. INTRODUCTION 



Miniaturization of fluidic devices from bench-top to palm-top size has been progressed 
considerably in recent years The advent of novel micro- and nano-manufacturing and 
fabrication techniques has equipped scientists and engineers the ability to manipulate the 
;ransport of micro/nano liters of fluid through the micro/nano channels in these devices 
The applications of such microfluidic devices are spread through fields like electronic- 
chip cooling, chemical synthesis, targeted cell isolation, bio-particle separation processes, 
chromatography, micr{>particle sorting, micro-reaction, micro- mixing, genomic/proteomic 



studies and others 
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Shrinking the size of table-top labs has led to the evolution 
of new kinds of 'on-chip' bio-assays such as lab-on-chip, blood-on-chip, cell-on-chip and 
neurons-on-chip [sl, [tI, l8| . The advantages of these microfluidic devices are manifold: they 
are portable, fast, affordable, accurate and energy efficient. The small sample and reagent 
volumes required and ease of use by non-experts also equip them to cater to 'point-of-care' 
needs. Fluid flow through channels of a few microns in size is a common feature of all these 
devices. 



The microchannels are, in many applications, coated with hydrophobic materials 



Hydrophobicity arises when the surface energy of the solid-liquid interface is 



usually measured from the contact angle made by a sessile drop on a surface 10|, lUj. Hy 



high and is 



drophobic coatings on the walls of the microchannels facilitate larger flow rates compared to 



hydrophilic counterparts for the same pressure drop as they offer less resistance to flow 12 1. 



Hydrophobic surfaces can also amplify electro-kinetic pumping, aid passive chaotic mixing 



and also mitigate the possibility of choking or adhesion of suspended analytes 13l-ll5[|. On 
hydrophobic surfaces, the traditional no-slip boundary condition is not valid and, instead, 
the fluid is modeled with a finite velocity at the wall. Navier proposed a generalized bound- 
ary condition to model the velocity of fiuid (u) tangential to the wall, by assuming it to be 



proportional to the surface shear stress 



16| 



dy 



where the slip-length /3 is the distance from the surface to the point where the linearly 
extrapolated velocity profile vanishes. The slip-length can be used to characterize the type 
of fiow in channels; if /3 = the fiow is stick-fiow (i.e. no slip), if /3 = 00 the fiow is 



plug-flow (e.g. shear free boundary) and any value of f3 between these two would represent 
partial slip-flow. In practical terms, the slip length at hydrophobic surfaces varies from a 
few nanometers to a few microns 17|. A slip length of up to 185 fim was reported in an 



experimental study [18 



which is comparable to the size of boundary layers in macroscopic 



regime. In the macroscopic regime the no-slip boundary condition {i.e. velocity of the 
:luid at solid surface is equal to the wall velocity) captures the physics of flow adequately 



9|, [lol. However, for flows in microchannels, partial-slip boundary conditions have to be 



applied at hydrophobic surfaces. In the remainder of this paper, we will use hydrophillic to 
refer to surfaces on which slip length is insignificant and the no-slip boundary condition is 
appropriate and hydrophobic for surfaces with a finite partial slip. 

Numerical simulation of flow through microchannels is important for understanding the 
underlying physics of these flows, as well as minimizing effort and expense of experiments, 
especially during the design and optimization of microfluidic devices. The range of numer- 
ical simulation methods spans continuum based computational fluid dynamics (CFD) to 
atomistic level molecular dynamics (MD). Modeling flow through microchannels based on 
the continuum assumption suffers from several drawbacks. First, for flows with high Knud- 



sen number (Kn) the continuum approximation may begin to fail [20j. Second, many of 
the flows include the presence of mesoscale particles such as DNA or individual biological 
cells. Treating the interaction of such second phase particles with the fluid medium and 
including the Brownian effects due to random thermal fluctuations becomes computation- 
ally prohibitive in CFD calculations. Microscopic modeling of above problems using MD 
is also not a practical choice as it is much too detailed and computationally expensive. As 
alternatives, discrete computational schemes like stochastic rotation dynamics, Brownian 
dynamics, lattice-Boltzmann method, smoothed particle dynamics, DPD have been devel- 



oped primarily for the spatio-temporal scales relevant to these situations j2l|. Among these 



discrete methods, DPD is a popular mesoscale scheme which bridges the gap between the 



macroscopic CFD and the microscopic MD 22 1. 



DPD was introduced to study the dynamics of complex fluids such as colloids, soft matter 



and polymers |23[. A cluster of atoms or molecules are considered to form a single particle 
in the DPD scheme. The positions of the particles are updated in a Lagrangian framework, 
and thus DPD can be viewed as a coarse-grained version of MD. In the initial version 
of the method, the DPD particles were treated as point masses and larger sized particles 



could be modeled only by binding several of the DPD particles together appropriately. A 



modified version of DPD was introduced later [2J] by treating particles as finite sized and 
solving for the concomitant rotational degrees of freedom. In this finite-size DPD (FDPD) 
model, non-central and rotational dissipative forces were considered in addition to the central 
conservative and non-conservative forces used in conventional DPD. 

The incorporation of appropriate wall boundary conditions for DPD has proved to be a 



challenge 
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26| . Discrete methods show density fiuctuations near the walls. While 



such density fiuctuations are realistic at the molecular level, they are considered spurious at 
the mesoscale continuum level at which DPD purports to model fiuids. An instantaneous 
wall boundary (IWB) model introduced by Ranjith et al. 27| proved useful for modeling 



Dations. Infiow and 



28|. 



no-slip and partial-slip with minimum near wall fiuid property perturl 
outfiow boundary conditions have also recently been modeled in DPD 

The hydrodynamics of the developing region inside channels with no-slip surfaces has been 
studied experimentally analytically 30 1 and numerically using CFD 31 1. Even though 
there have been considerable efforts to study the fluid transport in hydrophobic microchan- 
nels j^, [l3| in the fully develo ped reg ime, not much attention had been paid to the entrance 



effects with a few exceptions 



32 



Since the entrance length is proportional to Reynolds 



number Re, it can reasonably be ignored in low Re flows. While typical microflows are 



characterized by Re < 30 



34l |. in a few microfluidic applications 



micro-mixers etc the Re reaches the order of a few hundreds 
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ie micro heat-exchangers, 
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36|. Wall shear stress 



effects and velocity distributions vary signiflcantly at the entrance and these may eventually 



affect the separation efficiency of the microfluidic processes |32|. Moreover, the entrance 
region in hydrophobic channels is much longer than hydrophilic channels |33|. Reduction 
of entrance length of hydrophobic microchannels is very important for the design of some 
types of LOG devices. In the present work, we study the hydrodynamics of developing flow 
between two parallel plates with hydrophilic and hydrophobic surfaces using FDPD method 
for the flrst time. We also study the effect of hydrophilic patches in the entrance of the 
channel and their effect in reducing the developing length in microchannels. 

In section |lll analytical solutions for flow between two inflnite parallel plates with partial 
slip are summarized briefly following 30, Q, Is^]. In section UTTl the governing equations 
of the FDPD method are summarized. The details for the implementation of slip-wall 
boundary conditions and inflow and outflow boundary conditions are also presented. The 



FDPD simulation results are then compared with the analytical solution of the developing 
flow for both no-slip and slip flows in section IIVI The velocity proflles and the developing 
length compare well with the analytical solutions. We then study flow in a channel with 
hydrophobic walls with a short hydrophilic strip at inlet in section |Vl We flnd that the 
hydrophilic strip at the entrance shortens the development length signiflcantly. 



II. ANALYTICAL SOLUTION FROM CONTINUUM THEORY 



In this section, we briefly summarize the analytical solution of developing flow between 
parallel flat plates following Sparrow et al 30|, Chakraborty and Anand 33|, and Duan and 



Muzychka 37| closely. Consider pressure driven flow between two inflnitely long parallel 
plates separated hj H = 2h. The steady, incompressible flow is governed by the mass 
balance 

0, 



du 
dx 



dv 
dy 



and momentum balance 



du dv 
dx dy 



1 dp 
p dx 



V- 



d'^u 



(2) 



(3) 



equations in two dimensions, where u and v are velocities in x and y direction, p is the 



density and v is the kinematic viscosity of fluid. The linearized momentum equation 
of the form 



30j is 



d'^u 



V- 



u- 



_du 
dx 



V f du\ 



(4) 



dy '^\^y/y=h 
where u is the average cross sectional velocity. The analytical solution of the governing 



equation was obtained following Sparrow et al. [30| by assuming a no- slip boundary con- 
dition. Later, Chakraborty and Anand 33l|, and Duan and Muzychka 37| assumed that 
the fluid has a flnite velocity at the wall and is modeled by the Navier boundary condition 
(Eq. [1]) jssi Q . The slip- length /3 was modeled as a function of properties of the gas layer 
adjacent to the wall in case of flow between hydrophobic surface by 33| and as a function of 
Kn in case of fluid flow through microchannels for 0.001 < Kn < 0.1 by js^]. The governing 
equations are nondimensionalized using the hydraulic diameter Dh (which is 2H for parallel 
plates), half width h and average velocity u. Hence the dimensionless parameters are taken 
as ^ = x/(f>, T] = y/h, f3' = f3/h and U = u/u where = DyJ {puDh/ p). The dimensionless 
form of Navier boundary condition is given hj U = (3' {dU /drj). The dimensionless steady 



velocity profile is a function of both spatial coordinates ^ and rj and is given by 



U( c) - ^f^' 3(1 -r/^) >^ 2[a^cos(aj?7) - sin(aj)]exp(-16a^^0 
The eigenvalues aj satisfy 

The eigenvalues for partial slip 7^ 0) and no-slip = 0) cases are listed in the appendix. 



III. SIMULATION PROCEDURE AND BOUNDARY CONDITIONS 



A. Finite-size dissipative particle dynamics 



In this section, the formulation of FDPD model presented by Pan et al. 2^ is sum- 
marized. The domain of interest consists of DPD particles of finite size with a number 
density p. In this model, a set of molecules constitute an FDPD particle having a mass mi 
(i = 1 , . . . , A^) and mass moment of inertia /j (i = 1 , . . . , A^) . The degree of coarse-graining 
depends on the degree of spatio-temporal detail required. Each FDPD particle obeys New- 
ton's laws of motion and the translational motion is governed by the linear momentum 
equation 

where Vj and fj are, respectively, the velocity of and the force on the iih particle. As the 
particles are of finite size, the angular momentum equation is enforced by 

. dijJi 



'^XijVij X fij, (8) 



dt 

where (jJi is the angular velocity and fy is the effective force exerted on the ith particle 
by the neighboring jth particle, at a distance Tij = — r^- . The tangential forces are 
assumed to impart torques on the particles in proportion to the particle radii Ri and thus 
Xij = Ri/{Ri + Rj). The position, linear, and angular velocities of each particle is determined 
by the total force exerted by the surrounding particles within a certain finite cut-off radius 
Tc. In this scheme, the contribution from four types of forces (central (C), translational (T), 
rotational (R) and stochastic (5)) are considered 



The total force on the ith particle due to the surrounding particles is given by 

The total force on the ith particle may include external forces (f"^) from gravitational, 
magnetic or electro-osmotic forces, if any, f; = f; + f^. 

The central conservative repulsive force acting along the line connecting centers is taken 
to be 

fj = a,,r(r,,)e,„ (11) 

where aij is the repulsion parameter, = |rjj| and e^j = Vij/rij is a unit vector. An 
appropriate weight function T{rij) is selected such that the conservative force decreases 
monotonically to at r^j = Tc- Most DPD simulations employ the form of r(rjj) given by 

{1 _ lii if r < , , 

(12) 
if rij > Tc. 

The translational force is assumed to have central and non-central dissipative components 
given by 

= -7ii^r2(rij)(vij ■ eij)eij - -fij^T^{rij)[^rij - (v^j ■ eij)eij], (13) 

where 7jj'" and 'jij^ are the central and shear dissipation coefficients respectively. This 
frictional force attempts to reduce the relative velocity Vjj = Vj — between particles in 
both directions. The rotational dissipative force is taken to be of the form 

^5 = -7ii^r2(rij.)[ri^. X {\,jLJi + XjiUJj)]. (14) 

Finally, a stochastic force is also accounted through the expression given below 

fi^At = T{r,j)[a,j^tr[dW,j]^l + v^a./rfW^/] • e^,-, (15) 

where At is the time step, d = 2 for two-dimensional simulations and tr[dWij] is the trace 
of symmetric independent Wiener increment matrix dW ij while dWij"^ is its antisymmetric 
part. According to the fluctuation dissipation theorem, the random and dissipation coeffi- 



cients are related by cTjj'^ = ^j2kBT'^if and a if = a/2AJ^T7^. The stochastic forces and 
the dissipation forces together act to maintain a constant temperature during the simula- 
tions. 
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B. Wall boundary conditions 



As mentioned in the Introduction, DPD simulations have shown spurious density fluctu 



ations at walls when enforcing no-slip boundary conditions |25l. |26| . Recently, a new method 
of enforcing wall boundary conditions in FDPD simulations has shown substantial reduction 
in the density fluctuations |27|. In this method, when a fluid particle is within the range of 
influence of the wall, the particle interacts with the closest point on the wall as if there were 
a wall particle for that time step. The interaction of the wall particle and the fluid particle is 
separately specified. This proved to be a simple method to reduce spurious density variations 
as well as control slip at the wall. This method, referred to as the instantaneous wall par- 
ticle boundary (IWB), was shown to be a computationally efficient procedure for modeling 
impenetrable walls. The slip velocity at wall was tuned by controlling the lateral dissipative 
force component between fluid and wall along the direction tangential to the wall. This is 
achieved by changing the lateral dissipation coefficient = a(l — rp^/rc)^7pp'^, which acts 
only within a distance Tc from the wall. Here Vpu, is the distance between wall and fluid 
particle. As in the fluid particle-particle interactions, the dissipative and random coeffi- 



cients are related by ap^^ = ^J^lk^T^i^ ■ The slip-length increases as the slip modification 
parameter a is decreased. Thus boundary conditions ranging from no-slip to a large partial 
slip could be achieved by tuning a slip modification factor a described in that scheme {27! . 
We note that, some surfaces achieve superhydrophobicity by enhanced roughness which trap 
air pockets. In our simulations, we do not model roughness or the second phase fiuids at 
the wall. Instead we specify an effective slip velocity using the Navier condition given by 
Eq. ([H). 



C. Determination of slip-length of a hydrophobic surface 

A simulation box of 20rc x 30rc size taken to be periodic in the stream-wise direction and 
bound by IWB walls in y direction is used to estimate the slip-length /3. For a constant body 
force (f'^ = 0.01), the volume flow-rate per unit area Q is calculated. The /?' is estimated 



from the theoretical expression 



38 



Q 



slip 



Q 



no— slip J 



3/3' (16) 



where slip flow-rate QsUp for different hydrophobic surfaces were obtained by tuning the 
parameter a. As seen from Fig. [H the flow-rate increases with increasing /?'. Thus the 
flow-rate obtained for a certain applied body force increases with increasing partial slip in 



13, 



accordance with the experimental findings reported earlier 

The expression for fully developed velocity profile for slip flow in dimensionless form is 
given by [33|, 



6/3' + 2 6/3' + 2 

Furthermore, the value of /3' of a hydrophobic surface in the FDPD simulations is obtained 
by fitting the analytical velocity profile (Eq. [T7|) to the simulated fully developed velocity 
profile for slip flow, refer Fig. |2l The difference in /3' determined by both methods is less 
than 0.66% for the range of slip lengths considered. For a = 0.15 the slip-length is found 
to be /3 = 2.25rc and corresponding /3' = (2.25/10) = 0.225. Unless otherwise specified, this 
value is used to model all the hydrophobic surfaces mentioned in this work. The FDPD 
simulation of two long hydrophilic surfaces separated hj H = SOr^ and H = 20rc has been 



carried out to check the effect of the width on the flow-rate. The theoretical sca^ 



flow-rate per unit area {Q^nolsUpl Qno-sUp) for /3 = obtained from the expression 39| 



ing of 



/l2 / 


' _dv\ 




3^1 I 


^ dx J 


n 



(18) 



is 2.25 and that from FDPD simulation is 2.251. Moreover for slip flow {13 = 2.25rc) the 

)30rc 

' slip : 



flow-rate ratio from the simulation Qfilp /Ql^lp is 1.954 while theoretical prediction is 1.95. 



The FDPD scheme along with the IWB wall is thus able to capture the developed flow 
hydrodynamics of hydrophilic and hydrophobic parallel plates. 



D. Channel inflow and outflow conditions 



In order to compare the FDPD simulations with analytical results, the inflow has been 
ensured to be uniform at the inlet with a velocity ( n = 1) at unit temperature {ksT = 1). To 
maintain the number density of particles in the channel, the particles leaving the channel at 
the outflow at each time step are reintroduced at the inlet with random positions and random 
velocities. The angular and translational velocities are drawn from a uniform distribution 
in such a way that the system temperature (fc^T = 1) is not affected by the newly inserted 
particles. However, it is observed that the reintroduced particles experience forces only 
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from the fluid domain and decelerate. The inlet velocity profiles do not match the analytical 
results due to this deceleration. In order to overcome this problem, a set of fixed particles 
(with purely conservative interaction potential) are introduced at the inlet with the same 
number density of the fluid particles in the domain. These particles provide a balancing 
force for the particles being reintroduced to the channel and allow the inlet flow to be at 
the required uniform velocity. 

Similarly, the outflow boundary conditions require balancing forces from outside the fluid 
domain. In the absence of such balancing forces, the fluid accelerates near the outflow 
region. To mitigate this effect, particles are fixed outside the outlet of the fluid domain 
at the same number density as in the fluid domain. These downstream repulsive particles 
provide the requisite opposing force to the fluid particles at the outlet (schematically shown 
in Fig. |3]). The repulsive interaction forces exerted by these particles on the fluid particles 
in outlet is calculated by trial and error. The FDPD velocity profiles were found to match 
the analytical fully developed profiles to good accord for an inter-particle conservative force 
parameter value of a^o = napp with n = 1.2. 



IV. FDPD SIMULATIONS OF DEVELOPING FLOW 

We devote this section to study the ability of the FDPD model for simulating developing 
flows in channels with no-slip and partial slip boundary conditions. We consider a steady 
developing flow between two long parallel plates. The size of the 2D simulation domain 
is taken to be 20rc x AOOvc and filled with p = 3 particles per unit volume. The effective 
fluid viscosity is calculated to be = 1 ^^j. The short range interactions were calculated 
with a cut-off radius = 1. The Reynolds number Re = (puH/n) is calculated to be 
60 for the geometry and fluid properties under consideration. The wall and fluid particles 
are taken to be of the same size and thus the particle size coefficients are Ajj = Xji = |. 
The maximum repulsion parameter between fluid-fluid a^p = IbksT / p is chosen according 



to Groot and Warren 22[ and fluid-wall ap^ = 20 is taken Op' and 'ty' represent fluid 



and wall particles respectively) according to Ranjith et al. 27|. The central and shear 



coefficients of the dissipative and random forces are taken to be 7^ = 7^^ = 7^ = 4.5 and 
Gpp = app = Gp^ = 3. The time step was taken to be dt = 0.01. 



All particles in the domain are arranged randomly with zero initial velocity. A body 
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force is assumed to act on each particle, which accelerates through the domain. The 
net momentum (U) of domain increases from to a uniform value of 1. There domain is 
decomposed into 400 x 200 bins in x and y direction across the length and breadth of channel 
to obtain statistical averages of the velocity inside the domain. The component of velocity in 
the direction of flow was averaged over 2 x 10^ iterations to get statistically accurate results. 
The force is adjusted for each partial-slip boundary condition to maintain the flow rate 
of Q = 1. For the range of forces applied in the present simulations, a value of k = 1.2 (ratio 
of the interaction coefficient of the fluid particle and inflow and outflow boundary particles 
and inter-particle interaction coefficients) ensured that Q = 1 and the velocity profiles are 
very close to the analytical solution. It was observed that for k < 1.15 the particles close to 
the outlet accelerate and for k, > 1.25 they decelerate. In both cases the simulated velocity 
field did not match the analytical solution. 



A. Flow in a long hydrophilic channel 



A uniform velocity profile at inlet with average velocity u at the inlet transforms to a 
parabolic velocity profile at the outlet with a maximum velocity 1.5n. Within the developing 
region the velocity is a function of both x and y. When the fiow is fully developed, velocity 
profile is given by Eq. (fTTl) and remains same further downstream. 

The no-slip condition at the wall is obtained by modelin g th e solid boundary with a 



slip modification factor a = 3 as reported in our earlier work 27|. The infiow and outfiow 
boundary conditions were implemented as discussed in section IIIIDI The velocity profiles at 
different axial positions {C, = constant) are plotted in Fig. H] (a), along with the analytical 
solution. In Fig. 111(b), the velocity profiles at different heights (77 = constant) are extracted 
and compared with the analytical solution. The FDPD simulation and analytical solutions 
were found to be in good agreement. The channel length was chosen to be Lc = 20H to 
minimize end effects. Some end effects are apparent within H/A from the outlet in the form 
of velocity fluctuations (less than 12% of the maximum velocity). 
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B. Flow in a long hydrophobic channel 



The effective slip at the hydrophobic surface is obtained by choosing an appropriate 
parameter a as discussed in section IIII B[ Due to the non-zero velocity at the wall, the 
acceleration of the central laminar core is less compared to that of the no-slip case. The 
velocity gradients produced by partial-slip wall are smaller than the no-slip walls. Due 
to this, the developing length Lg for superhydrophobic microchannels is greater than in 



hydrophilic channels 



33|. 



The velocity profiles obtained through the FDPD simulation is compared to the ana- 
lytical solution given by Eq. (j5]). The theoretical and computational results are in good 
agreement as can be seen in Fig. [5l This simulation shows that the FDPD scheme, in 
combination with IWB wall, can capture the hydrodynamics of the steady non-uniform de- 
veloping region of partial-slip flow accurately. The effect of a no-slip region at the entrance 
on the hydrodynamics of fluid flow between two hydrophobic surfaces is discussed in the 
next section. 



V. ENTRANCE REGION WITH A HYDROPHILIC STRIP 

The developing length or entrance length is the stream-wise distance from the inlet to 
the point at which the boundary layers formed on both the walls merge. In an engineering 
sense, this is quantified by the distance along the flow direction at which the centerline 
velocity reaches 99% of the maximum fully developed velocity. There are several empirical 
correlations for the entrance length as a function of the Reynolds number for no-slip channels 
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371 ]. For moderate Reynolds numbers, the developing length constitutes a considerable 
portion of a miniaturized LOG device, and hence a reduction in the entrance length of the 
microchannel is highly desirable. For partial-slip flow, the development of the flow field 
occurs over a much greater length compared to no-slip wall boundary conditions. This is 
because the acceleration of the central core for a no-slip wall is greater than that of a slip- 
wall for the same Re. The greater shear stress at the wall for a hydrophilic surface enables 
the presence of the wall to be felt inside the domain over a shorter distance compared to 
their superhydrophobic counterparts. The latter have a smaller magnitude of the velocity 
gradient near the walls. We explore the hydrodynamics of flow between two long parallel 

12 



superhydrophobic surfaces with a hydrophilic strip at the inlet in this section. 

Experiments on the entrance hydrodynamics in microchannels with no-shp boundary 
conditions with aspect ratio H/W = 1 {W being the width of the channel) were recently 



reported by Ahmad and Hassan [29|] for a hydraulic radius {Dh = H) ranging from lOO/xm 
to 500;um over a range of Re numbers from 0.5 to 200. The entrance length of a hydrophilic 
microchannel with Dh = 200/im and Re = 60 is interpolated from the empirical correlation 



obtained from their experimental data (Eg. (6) of Ref. 29|]) is 943 fim. Using an analyt 



ical approach, Chakraborty and Anand 



33| have presented a correlation that relates the 



developing length and Reynolds number of slip surfaces, 

0.044i?e(l + 1.675/3' + 2.3125/3'^). (19) 



Le 0.63 . nnAAn : , r^^ o> : o 01 oc: o/2n 



Dh 0.035i?e 

Here, Dh = H, for a given slip length (/?'). However, to the best of our knowledge there is no 
experimental data available till date on the hydrodynamics of entrance region of microchan- 
nels with hydrophobic surfaces. Hence, the above analytical solution in EqnJT9] is used to 
determine the developing length under partial-slip conditions. The ratio of the hydrody- 
namic development length for hydrophobic (partial-slip) to hydrophilic (no-slip) surfaces is 
calculated to be 

■ T 13' =0.225 ' 

1.46. 



r/9'=0 , 

Thus the developing lengths of hydrophobic channel is about 50% larger than hydrophilic 
channel for the same Re. It was reported in Ref. [l^ that, a nano-turf created by coating 
a surface with Teflon can produce a hydrophobic surface with /3 ^ 20 fim. For such a 
hydrophobic material with /?' = 0.225 (/3 = 22.5 fim, h = Dh/2 = 100 fim), and a typical 
Re = 60 would have a Lg ^ 1.46 x 943 fim = 1377 fim. For an 'on-chip' device which is 
already small in size, this undesirable entrance region would constitute a major portion of 
the channel. This effect increases as the number of parallel microchannels (which is typical 
for most LOG devices) increases to achieve large throughput. 

In this section the effect of introduction of a no-slip strip having length at the inlet, just 
before the hydrophobic surface (schematically shown in Fig. [6]) is discussed. The velocity 
profiles obtained through the DPD simulations closely follows the analytical solution, so Eq. 
(fT9|) was used to estimate Lg for slip flow. Thus for a hydrophobic surface with effective slip 
length /3' = 0.225 and Re = 60 has a developing length of SSr^. Fixing a hydrophilic strip 
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TABLE I. The entrance length for mixed hydrophobic-hydrophihc channels with different 
lengths of hydrophilic strips P at the inlet. 



SI. No. 






1 


Or, 


83rc 


2 


5rc 


50re 


3 


lOre 


30rc 


4 


15rc 


55rc 


5 


20rc 


65rc 


6 


25re 


lOOrc 



accelerates the central core faster compared with hydrophobic wall as large shear gradients 
are formed near no-slip surfaces. 

To study the influence of the no-slip strips, a number of DPD simulations were carried 
out by varying the initial hydrophilic strip length from Sr, to 25rc at Re = 60. As the 
stream-wise velocity gradient is maximum near the wall, the entrance length of a mixed 
hydrophilic-hydrophobic surface was found by monitoring the velocity near the wall. The 
transition between the hydrophilic and hydrophobic surface results in velocity fluctuations. 
These fluctuations are minimum along the centerline (77 = 0) and greatest at the wall. 
So the developing length with a hydrophilic strip (Lg) at the entrance, is estimated by 
determining the distance from the inlet to the point where the axial velocity for 77 = 0.9 
becomes constant. The effect of the length of the inlet hydrophilic strip on the developing 
length is given in TablelH It was found that for a lOrc hydrophilic strip, the developing length 
reduces from 83rc to SOvc (see Fig. [7]) and the percentage reduction in the developing length is 
(^^T~^) ~ 66%, although the portion of hydrophobic surface replaced by hydrophihc surface 
is only {l^ / Le) ~ | ~ 12.5%. Thus, the combination of hydrophilic-hydrophobic surfaces 
drastically reduced the developing length Lf,. The full development of velocity profile of such 
an arrangement takes place over a shorter distance than that of pure hydrophobic surfaces, 
as shown in Fig. [HI This finding is expected to be beneficial for the design of microfluidic 
devices and to optimize the size of LOG devices with hydrophobic channels. It is noteworthy 
that, from the analytical solution (Eq. [5]) the central line velocity at the end of hydrophilic 
strip (x = lOrc) is 1.244 and the fully developed centerline velocity of hydrophobic surface is 
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1.3. We infer that, the should be long enough to accelerate the central core of hydrophilic 
region nearer to the developed velocities of the hydrophobic surfaces. If the is too long 
the center region accelerates, and take a longer distance to decelerate and reach a developed 
state as shown in Fig. MJo). Conversely, if the l"^ is short the acceleration of central core is 
not enough to reach the uniform fully developed state and may take a longer distance, as 
seen in Fig. [8](e). 

VI. SUMMARY AND CONCLUSIONS 

In this work, a modified DPD method has been shown to effectively capture the hy- 
drodynamics of developing flows in microchannels with no-slip and partial-slip boundaries. 
The simulations with no-slip and partial slip wall boundary conditions were shown to have 
excellent agreement with analytical results. A new method to model inflow boundary con- 
dition is proposed to obtain a uniform inlet velocity profile. Similarly, the outflow boundary 
conditions were modified to prevent the fluid from accelerating out of the domain. The 
presence of a small hydrophilic strip at the inlet of a hydrophobic microchannel was found 
to significantly reduce the development length of the remaining hydrophobic channel. This 
finding can potentially be used by the designers of LOG devices to optimize the size of the 
microfluidic devices. 
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APPENDIX 



TABLE II. The eigenvalues obtained from Eq. [6] used to calculate the analytical solution of velocity. 



i 


13=0 


/3'=0.225 


1 


4.4934 


3.8666 


2 


7.7253 


6.8198 


3 


10.9041 


9.8327 


4 


14.0662 


12.8903 


5 


17.2208 


15.9749 


6 


20.3713 


19.0758 


7 


23.5195 


22.1871 


8 


26.6661 


25.3054 


9 


29.8116 


28.4286 


10 


32.9564 


31.5552 


11 


36.1006 


34.6845 


12 


39.2444 


37.8157 


13 


42.3879 


40.9485 


14 


45.5311 


44.0826 


15 


48.6741 


47.2176 


16 


51.8170 


50.3534 


17 


54.9597 


53.4898 


18 


58.1023 


56.6269 


19 


61.2447 


59.7644 


20 


64.3871 


62.9023 


21 


67.5294 


66.0406 


22 


70.6717 


69.1791 


23 


73.8139 


72.3180 


24 


76.9560 


75.4570 


25 


80.0981 


78.5963 
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FIG. 1. Variation of the flow-rate Q with slip-length f3 for a constant pressure gradient. 




FIG. 2. The fitted analytical solution following 
velocity profiles (shown with markers). 



33l | (shown with lines) over FDPD simulated 
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FIG. 3. Schematic representation of the inflow and outflow boundaries. 




FIG. 4. The velocity profiles in a hydrophilic microchannel = 0) along (a) span-wise = 
constant) and (b) stream-wise {rj = constant) directions. 
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FIG. 5. The velocity profiles in a hydrophobic microchannel (/?' = 0.225) along (a) = constant) 
and (b) stream-wise {r] = constant) directions. 
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FIG. 6. Schematic sketch of a hydrophobic channel with a hydrophilic strip of length at the 
inlet. 
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FIG. 7. Comparison of the simulated velocity profile of hydrophobic channel with hydrophilic 
inlet strip of length P = lOrc with the fully developed theoretical velocity profile at rj = 0.9 from 
Eq. (HH). 
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FIG. 8. Comparison of the simulated velocity profiles (shown with markers) at different locations 
from the inlet in a mixed hydrophobic channel for various lengths of hydrophilic strips at the inlet: 
(a) = Ore, (b) I' = 20rc, (c) = ISr^, (d) I' = lOr^, (e) = br^. The profiles at x = are 
marked A, x = lOrc by B, x = 30rc by C, x = SSr^ by D and x = 65rc by E in the figures. The 
fully developed analytical velocity profile for /3 = 2.25rc is shown in each case using a solid line for 



comparison. 



